! modal_aero_wateruptake.F90


!----------------------------------------------------------------------
!BOP
!
! !MODULE: modal_aero_wateruptake --- modal aerosol mode merging (renaming)
!
! !INTERFACE:
   module modal_aero_wateruptake

! !USES:
   use shr_kind_mod, only : r8 => shr_kind_r8
   use cam_logfile,  only: iulog

   implicit none
   private
   save
                                                                                                                             
!  !PUBLIC MEMBER FUNCTIONS:
   public modal_aero_wateruptake_sub
                                                                                                                             
! !PUBLIC DATA MEMBERS:
!  currently none

                                                                                                                             
! !DESCRIPTION: This module implements ...
!
! !REVISION HISTORY:
!
!   RCE 07.04.13:  Adapted from MIRAGE2 code
!
!EOP
!----------------------------------------------------------------------
!BOC
! list private module data here
!EOC
!----------------------------------------------------------------------
                                                                                                                             
                                                                                                                             
   contains
                                                                                                                             
                                                                                                                             
!----------------------------------------------------------------------
      subroutine modal_aero_wateruptake_sub(                &
                 lchnk, ncol, nstep,                        &
                 iwaterup_flag, loffset,                    &
                 aero_mmr_flag, h2o_mmr_flag,               &
                 deltat, h2ommr, t, pmid, pdel, cldn,       &
                 raer, raertend, dotend, qaerwat,           &
                 dgncur_a, dgncur_awet, wetdens             &
#ifdef OBSRH
                 , obsrh                                    &
#endif
                                                            )
! following are local for now -- wait and see
!                maer, naer, wetrad, density                )

!-----------------------------------------------------------------------
!
! Purpose: Compute aerosol wet radius
!
! Method:  Kohler theory
!
! Author:  S. Ghan
!
!-----------------------------------------------------------------------
      use modal_aero_data

      use mo_constants,  only: pi
      use pmgrid,        only: plat, plon
      use ppgrid,        only: begchunk, endchunk, pcols, pver
      use physconst,     only: cpair, epsilo, gravit, mwdry, mwh2o,   &
                               rair, rga, rhoh2o, rh2o, latvap
      use constituents,  only: pcnst, qmin
      use wv_saturation, only: aqsat, qsat_water
      use phys_grid,     only: get_lat_all_p, get_lon_all_p

      use abortutils,    only : endrun
      use cam_history,   only : outfld
      use spmd_utils,    only : masterproc

      implicit none

      integer,  intent(in)  :: lchnk              ! chunk index
      integer,  intent(in)  :: ncol               ! number of columns
      integer,  intent(in)  :: nstep              ! time step
      integer,  intent(in)  :: iwaterup_flag       
                               ! identifies call location from typhysbc (1,2)
      integer,  intent(in)  :: loffset            ! offset applied to modal aero "pointers"

      logical,  intent(in)  :: aero_mmr_flag      ! if .true.,  aerosol q are kg/kg-air
                                                  ! if .false., aerosol q are mol/mol-air
      logical,  intent(in)  :: h2o_mmr_flag       ! if .true.,  h2ommr is kg/kg-air
      logical,  intent(inout)::dotend(pcnst)
                               ! identifies species for which tendencies are computed

      real(r8), intent(in)  :: deltat             ! time step (s)
      real(r8), intent(in)  :: h2ommr(pcols,pver) ! layer specific humidity
      real(r8), intent(in)  :: t(pcols,pver)      ! layer temperatures (K)
      real(r8), intent(in)  :: pmid(pcols,pver)   ! layer pressure (Pa)
      real(r8), intent(in)  :: pdel(pcols,pver)   ! layer pressure thickness (Pa)
      real(r8), intent(in)  :: cldn(pcols,pver)   ! layer cloud fraction (0-1)
      real(r8), intent(in)  :: raer(pcols,pver,pcnst)
                               ! aerosol species MRs (kg/kg and #/kg)
      real(r8), intent(inout)::raertend(pcols,pver,pcnst)
                               ! aerosol MR tendencies (kg/kg/s)
                               ! only defined for aerosol water
      real(r8), intent(out)   :: qaerwat(pcols,pver,ntot_amode)
      real(r8), intent(in)    :: dgncur_a(pcols,pver,ntot_amode)
      real(r8), intent(out)   :: dgncur_awet(pcols,pver,ntot_amode)
      real(r8), intent(out)   :: wetdens(pcols,pver,ntot_amode)

! following are local for now -- wait and see
!     real(r8), intent(out) :: maer(pcols,pver,ntot_amode)
!                              ! aerosol wet mass MR (including water) (kg/kg-air)
!     real(r8), intent(out) :: naer(pcols,pver,ntot_amode)
!                              ! aerosol number MR (bounded!) (#/kg-air)
!     real(r8), intent(out) :: wetrad(pcols,pver,ntot_amode)  
!                              ! wet radius of aerosol (m)
!     real(r8), intent(out) :: density(pcols,pver,ntot_amode)
!                              ! wet mean density of aerosol (kg/m3)

!     local variables

      integer i,k,m
      integer icol_diag
      integer l ! species index
      integer lmass ! pointer for aerosol mass
      integer ltype ! pointer for aerosol type
!     integer lwater ! pointer for water on aerosol
      integer  lat(pcols), lon(pcols)      ! lat,lon indices

      real(r8) density_water                   ! density of water (kg/m3)
      real(r8) drydens(ntot_amode)   ! dry particle density  (kg/m^3)
      real(r8) drymass(ntot_amode)   ! single-particle-mean dry mass  (kg)
      real(r8) dryrad(pcols,pver,ntot_amode) ! dry volume mean radius of aerosol (m)
      real(r8) dryvol(ntot_amode)    ! single-particle-mean dry volume (m3)
      real(r8) dryvolmr(ntot_amode)  ! volume MR for aerosol mode (m3/kg)
      real(r8) duma, dumb
      real(r8) es(pcols,pver)        ! saturation vapor pressure (Pa)
      real(r8) hygro(ntot_amode)     ! volume-weighted mean hygroscopicity (--)
      real(r8) hystfac(ntot_amode)   ! working variable for hysteresis
      real(r8) pi43
      real(r8) qs(pcols,pver)        ! saturation specific humidity
      real(r8) qwater                ! aerosol water MR
      real(r8) rh(pcols,pver)        ! relative humidity (0-1)
      real(r8) third
      real(r8) v2ncur_a(pcols,pver,ntot_amode)
      real(r8) wtrvol(ntot_amode)    ! single-particle-mean water volume in wet aerosol (m3)
      real(r8) wetvol(ntot_amode)    ! single-particle-mean wet volume (m3)

      real(r8) :: maer(pcols,pver,ntot_amode)
                              ! aerosol wet mass MR (including water) (kg/kg-air)
      real(r8) :: naer(pcols,pver,ntot_amode)
                              ! aerosol number MR (bounded!) (#/kg-air)
      real(r8) :: wetrad(pcols,pver,ntot_amode)  
                              ! wet radius of aerosol (m)

      character(len=3) :: trnum       ! used to hold mode number (as characters)

!-----------------------------------------------------------------------



! set the dotend's
! currently, this is one of several routines called from aerosol_wet_intr
! so DO NOT re-initialize dotend=.false. and raertend=0.0
!     dotend(:) = .false.
      do m=1,ntot_amode
!        lwater = lwaterptr_amode(m) - loffset
!        dotend(lwater) = .true.
!        raertend(1:ncol,:,lwater) = 0.0
      end do


      third=1./3.
      pi43 = pi*4.0/3.0
      density_water = rhoh2o   ! is (kg/m3)

! compute size-related factors
!    when numptr_amode(m)=0, earlier code set dryrad = prescribed volume.
!    this is no longer needed because
!           dryrad(i,k,m) = (1.0/v2ncur_a(i,k,m)/pi43)**third
!    works in all cases
      do m=1,ntot_amode
           hystfac(m) = 1.0 / max( 1.0e-5_r8,   &
                              (rhdeliques_amode(m)-rhcrystal_amode(m)) )
      enddo

! main loops over i, k

!     call aqsat( t, pmid, es, qs, pcols, ncol, pver, 1, pver )
      do k=1,pver
      do i=1,ncol

         qs(i,k)=qsat_water(t(i,k),pmid(i,k))

         if ( h2o_mmr_flag ) then
            rh(i,k) = h2ommr(i,k)/qs(i,k)
         else
            rh(i,k) = h2ommr(i,k)*mwh2o/(mwdry*qs(i,k))
         end if
         rh(i,k) = max(rh(i,k),0.0_r8)
         rh(i,k) = min(rh(i,k),0.98_r8)
         if (cldn(i,k) .lt. 1.0_r8) then
           rh(i,k) = (rh(i,k) - cldn(i,k)) / (1.0_r8 - cldn(i,k))  ! clear portion
         end if
         rh(i,k) = max(rh(i,k),0.0_r8)
           

!     compute dryvolmr, maer, naer for each mode
         do m=1,ntot_amode

            maer(i,k,m)=0.
            dryvolmr(m)=0.
            hygro(m)=0.
            do l = 1, nspec_amode(m)
               lmass = lmassptr_amode(l,m) - loffset
               ltype = lspectype_amode(l,m)
               if ( aero_mmr_flag ) then
                  duma = raer(i,k,lmass)
               else
                  duma = raer(i,k,lmass)*(specmw_amode(ltype)/mwdry)
               end if
               maer(i,k,m) = maer(i,k,m) + duma
               dumb = duma/specdens_amode(ltype)
               dryvolmr(m) = dryvolmr(m) + dumb
               hygro(m) = hygro(m) + dumb*spechygro(ltype)
            enddo
            if (dryvolmr(m) > 1.0e-30_r8) then
               hygro(m) = hygro(m)/dryvolmr(m)
            else
               hygro(m) = spechygro( lspectype_amode(1,m) )
            end if

!     naer = aerosol number (#/kg)
!     the new v2ncur_a replaces old coding here
            v2ncur_a(i,k,m) = 1. / ( (pi/6.)*                            &
                (dgncur_a(i,k,m)**3.)*exp(4.5*alnsg_amode(m)**2.) )
            naer(i,k,m) = dryvolmr(m)*v2ncur_a(i,k,m)
         enddo   !m=1,ntot_amode

!     compute mean (1 particle) dry volume and mass for each mode
!     old coding is replaced because the new (1/v2ncur_a) is equal to
!        the mean particle volume
!     also moletomass forces maer >= 1.0e-30, so (maer/dryvolmr)
!        should never cause problems (but check for maer < 1.0e-31 anyway)
         do m=1,ntot_amode
            if (maer(i,k,m) .gt. 1.0e-31) then
               drydens(m) = maer(i,k,m)/dryvolmr(m)
            else
               drydens(m) = 1.0
            end if
            dryvol(m) = 1.0/v2ncur_a(i,k,m)
            drymass(m) = drydens(m)*dryvol(m)
            dryrad(i,k,m) = (dryvol(m)/pi43)**third
         enddo

!     compute wet radius for each mode
         do m=1,ntot_amode
            call modal_aero_kohler(                          &
                    dryrad(i,k,m), hygro(m), rh(i,k),        &
                    wetrad(i,k,m), 1, 1                      )

            wetrad(i,k,m)=max(wetrad(i,k,m),dryrad(i,k,m))
            dgncur_awet(i,k,m) = dgncur_a(i,k,m)*   &
                                     (wetrad(i,k,m)/dryrad(i,k,m))
            wetvol(m) = pi43*wetrad(i,k,m)*wetrad(i,k,m)*wetrad(i,k,m)
            wetvol(m) = max(wetvol(m),dryvol(m))
            wtrvol(m) = wetvol(m) - dryvol(m)
            wtrvol(m) = max( wtrvol(m), 0.0_r8 )

!     apply simple treatment of deliquesence/crystallization hysteresis
!     for rhcrystal < rh < rhdeliques, aerosol water is a fraction of
!     the "upper curve" value, and the fraction is a linear function of rh
            if (rh(i,k) < rhcrystal_amode(m)) then
               wetrad(i,k,m) = dryrad(i,k,m)
               wetvol(m) = dryvol(m)
               wtrvol(m) = 0.0_r8
            else if (rh(i,k) < rhdeliques_amode(m)) then
               wtrvol(m) = wtrvol(m)*hystfac(m)   &
                                    *(rh(i,k) - rhcrystal_amode(m))
               wtrvol(m) = max( wtrvol(m), 0.0_r8 )
               wetvol(m) = dryvol(m) + wtrvol(m)
               wetrad(i,k,m) = (wetvol(m)/pi43)**third
            end if

!     compute aer. water tendency = (new_water - old_water)/deltat
!     [ either (kg-h2o/kg-air/s) or (mol-h2o/mol-air/s) ]
!           lwater = lwaterptr_amode(m) - loffset
            if ( aero_mmr_flag ) then
               duma = 1.0_r8
            else
#if (defined MIMIC_CAM3)
               duma = mwdry/18.0_r8
#else
               duma = mwdry/mwh2o
#endif
            end if
            qwater = density_water*naer(i,k,m)*wtrvol(m)*duma

!     old_water (after modal_aero_calcsize) is 
!           qwater_old = raer(i,k,lwater) + raertend(i,k,lwater)*deltat
!     and water tendency is
!           raertend(i,k,lwater) = raertend(i,k,lwater)   &
!                                + (qwater - qwater_old)/deltat
!     which is equivalent to
!           raertend(i,k,lwater) = (qwater - raer(i,k,lwater))/deltat
            qaerwat(i,k,m) = qwater

!     compute aerosol wet density (kg/m3)
            if (wetvol(m) > 1.0e-30_r8) then
               wetdens(i,k,m) = (drymass(m) + density_water*wtrvol(m))/wetvol(m)
            else
               wetdens(i,k,m) = specdens_amode( lspectype_amode(1,m) )
            end if

         enddo


      end do   ! "i=1,ncol"
      end do   ! "k=1,pver"


! output to history
      do m = 1, ntot_amode
         write( trnum, '(i3.3)' ) m
         call outfld( 'wat_a'//trnum(3:3),  qaerwat(:,:,m),     pcols, lchnk)
! note - eventually we should change these from "dgnd_a0N" to "dgnd_aN"
         call outfld( 'dgnd_a'//trnum(2:3), dgncur_a(:,:,m),    pcols, lchnk)
         call outfld( 'dgnw_a'//trnum(2:3), dgncur_awet(:,:,m), pcols, lchnk)
      end do


      return
      end subroutine modal_aero_wateruptake_sub


!-----------------------------------------------------------------------
      subroutine modal_aero_kohler(   &
          rdry_in, hygro, s, rwet_out, im, imx )

! calculates equlibrium radius r of haze droplets as function of
! dry particle mass and relative humidity s using kohler solution
! given in pruppacher and klett (eqn 6-35)

! for multiple aerosol types, assumes an internal mixture of aerosols

      implicit none

! arguments
      integer :: im         ! number of grid points to be processed
      integer :: imx        ! dimensioned number of grid points
      real(r8) :: rdry_in(imx)    ! aerosol dry radius (m)
      real(r8) :: hygro(imx)      ! aerosol volume-mean hygroscopicity (--)
      real(r8) :: s(imx)          ! relative humidity (1 = saturated)
      real(r8) :: rwet_out(imx)   ! aerosol wet radius (m)

! local variables
      integer, parameter :: imax=200
      integer :: i, n, nsol

      real(r8) :: a, b
      real(r8) :: p40(imax),p41(imax),p42(imax),p43(imax) ! coefficients of polynomial
      real(r8) :: p30(imax),p31(imax),p32(imax) ! coefficients of polynomial
      real(r8) :: p
      real(r8) :: r3, r4
      real(r8) :: r(imx)        ! wet radius (microns)
      real(r8) :: rdry(imax)    ! radius of dry particle (microns)
      real(r8) :: ss            ! relative humidity (1 = saturated)
      real(r8) :: slog(imax)    ! log relative humidity
      real(r8) :: vol(imax)     ! total volume of particle (microns**3)
      real(r8) :: xi, xr

      complex(r8) :: cx4(4,imax),cx3(3,imax)

      real(r8), parameter :: eps = 1.e-4
      real(r8), parameter :: mw = 18.
      real(r8), parameter :: pi = 3.14159
      real(r8), parameter :: rhow = 1.
      real(r8), parameter :: surften = 76.
      real(r8), parameter :: tair = 273.
      real(r8), parameter :: third = 1./3.
      real(r8), parameter :: ugascon = 8.3e7


!     effect of organics on surface tension is neglected
      a=2.e4*mw*surften/(ugascon*tair*rhow)

      do i=1,im
           rdry(i) = rdry_in(i)*1.0e6   ! convert (m) to (microns)
           vol(i) = rdry(i)**3          ! vol is r**3, not volume
           b = vol(i)*hygro(i)

!          quartic
           ss=min(s(i),1.-eps)
           ss=max(ss,1.e-10_r8)
           slog(i)=log(ss)
           p43(i)=-a/slog(i)
           p42(i)=0.
           p41(i)=b/slog(i)-vol(i)
           p40(i)=a*vol(i)/slog(i)
!          cubic for rh=1
           p32(i)=0.
           p31(i)=-b/a
           p30(i)=-vol(i)
      end do


       do 100 i=1,im

!       if(vol(i).le.1.e-20)then
        if(vol(i).le.1.e-12)then
           r(i)=rdry(i)
           go to 100
        endif

        p=abs(p31(i))/(rdry(i)*rdry(i))
        if(p.lt.eps)then
!          approximate solution for small particles
           r(i)=rdry(i)*(1.+p*third/(1.-slog(i)*rdry(i)/a))
        else
           call makoh_quartic(cx4(1,i),p43(i),p42(i),p41(i),p40(i),1)
!          find smallest real(r8) solution
           r(i)=1000.*rdry(i)
           nsol=0
           do n=1,4
              xr=real(cx4(n,i))
              xi=imag(cx4(n,i))
              if(abs(xi).gt.abs(xr)*eps) cycle  
              if(xr.gt.r(i)) cycle  
              if(xr.lt.rdry(i)*(1.-eps)) cycle  
              if(xr.ne.xr) cycle  
              r(i)=xr
              nsol=n
           end do  
           if(nsol.eq.0)then
              write(iulog,*)   &
               'ccm kohlerc - no real(r8) solution found (quartic)'
              write(iulog,*)'roots =', (cx4(n,i),n=1,4)
              write(iulog,*)'p0-p3 =', p40(i), p41(i), p42(i), p43(i)
              write(iulog,*)'rh=',s(i)
              write(iulog,*)'setting radius to dry radius=',rdry(i)
              r(i)=rdry(i)
!             stop
           endif
        endif

        if(s(i).gt.1.-eps)then
!          save quartic solution at s=1-eps
           r4=r(i)
!          cubic for rh=1
           p=abs(p31(i))/(rdry(i)*rdry(i))
           if(p.lt.eps)then
              r(i)=rdry(i)*(1.+p*third)
           else
              call makoh_cubic(cx3,p32,p31,p30,im)
!             find smallest real(r8) solution
              r(i)=1000.*rdry(i)
              nsol=0
              do n=1,3
                 xr=real(cx3(n,i))
                 xi=imag(cx3(n,i))
                 if(abs(xi).gt.abs(xr)*eps) cycle  
                 if(xr.gt.r(i)) cycle  
                 if(xr.lt.rdry(i)*(1.-eps)) cycle  
                 if(xr.ne.xr) cycle  
                 r(i)=xr
                 nsol=n
              end do  
              if(nsol.eq.0)then
                 write(iulog,*)   &
                  'ccm kohlerc - no real(r8) solution found (cubic)'
                 write(iulog,*)'roots =', (cx3(n,i),n=1,3)
                 write(iulog,*)'p0-p2 =', p30(i), p31(i), p32(i)
                 write(iulog,*)'rh=',s(i)
                 write(iulog,*)'setting radius to dry radius=',rdry(i)
                 r(i)=rdry(i)
!                stop
              endif
           endif
           r3=r(i)
!          now interpolate between quartic, cubic solutions
           r(i)=(r4*(1.-s(i))+r3*(s(i)-1.+eps))/eps
        endif

  100 continue

! bound and convert from microns to m
      do i=1,im
         r(i) = min(r(i),30._r8) ! upper bound based on 1 day lifetime
         rwet_out(i) = r(i)*1.e-6
      end do

      return
      end subroutine modal_aero_kohler


!-----------------------------------------------------------------------
      subroutine makoh_cubic( cx, p2, p1, p0, im )
!
!     solves  x**3 + p2 x**2 + p1 x + p0 = 0
!     where p0, p1, p2 are real
!
      integer, parameter :: imx=200
      integer :: im
      real(r8) :: p0(imx), p1(imx), p2(imx)
      complex(r8) :: cx(3,imx)

      integer :: i
      real(r8) :: eps, q(imx), r(imx), sqrt3, third
      complex(r8) :: ci, cq, crad(imx), cw, cwsq, cy(imx), cz(imx)

      save eps
      data eps/1.e-20/

      third=1./3.
      ci=dcmplx(0.,1.)
      sqrt3=sqrt(3.)
      cw=0.5*(-1+ci*sqrt3)
      cwsq=0.5*(-1-ci*sqrt3)

      do i=1,im
      if(p1(i).eq.0.)then
!        completely insoluble particle
         cx(1,i)=(-p0(i))**third
         cx(2,i)=cx(1,i)
         cx(3,i)=cx(1,i)
      else
         q(i)=p1(i)/3.
         r(i)=p0(i)/2.
         crad(i)=r(i)*r(i)+q(i)*q(i)*q(i)
         crad(i)=sqrt(crad(i))

         cy(i)=r(i)-crad(i)
         if (abs(cy(i)).gt.eps) cy(i)=cy(i)**third
         cq=q(i)
         cz(i)=-cq/cy(i)

         cx(1,i)=-cy(i)-cz(i)
         cx(2,i)=-cw*cy(i)-cwsq*cz(i)
         cx(3,i)=-cwsq*cy(i)-cw*cz(i)
      endif
      enddo

      return
      end subroutine makoh_cubic


!-----------------------------------------------------------------------
      subroutine makoh_quartic( cx, p3, p2, p1, p0, im )

!     solves x**4 + p3 x**3 + p2 x**2 + p1 x + p0 = 0
!     where p0, p1, p2, p3 are real
!
      integer, parameter :: imx=200
      integer :: im
      real(r8) :: p0(imx), p1(imx), p2(imx), p3(imx)
      complex(r8) :: cx(4,imx)

      integer :: i
      real(r8) :: third, q(imx), r(imx)
      complex(r8) :: cb(imx), cb0(imx), cb1(imx),   &
                     crad(imx), cy(imx), czero


      czero=cmplx(0.0,0.0)
      third=1./3.

      do 10 i=1,im

      q(i)=-p2(i)*p2(i)/36.+(p3(i)*p1(i)-4*p0(i))/12.
      r(i)=-(p2(i)/6)**3+p2(i)*(p3(i)*p1(i)-4*p0(i))/48.   &
       +(4*p0(i)*p2(i)-p0(i)*p3(i)*p3(i)-p1(i)*p1(i))/16

      crad(i)=r(i)*r(i)+q(i)*q(i)*q(i)
      crad(i)=sqrt(crad(i))

      cb(i)=r(i)-crad(i)
      if(cb(i).eq.czero)then
!        insoluble particle
         cx(1,i)=(-p1(i))**third
         cx(2,i)=cx(1,i)
         cx(3,i)=cx(1,i)
         cx(4,i)=cx(1,i)
      else
         cb(i)=cb(i)**third

         cy(i)=-cb(i)+q(i)/cb(i)+p2(i)/6

         cb0(i)=sqrt(cy(i)*cy(i)-p0(i))
         cb1(i)=(p3(i)*cy(i)-p1(i))/(2*cb0(i))

         cb(i)=p3(i)/2+cb1(i)
         crad(i)=cb(i)*cb(i)-4*(cy(i)+cb0(i))
         crad(i)=sqrt(crad(i))
         cx(1,i)=(-cb(i)+crad(i))/2.
         cx(2,i)=(-cb(i)-crad(i))/2.

         cb(i)=p3(i)/2-cb1(i)
         crad(i)=cb(i)*cb(i)-4*(cy(i)-cb0(i))
         crad(i)=sqrt(crad(i))
         cx(3,i)=(-cb(i)+crad(i))/2.
         cx(4,i)=(-cb(i)-crad(i))/2.
      endif
   10 continue

      return
      end subroutine makoh_quartic

!----------------------------------------------------------------------

   end module modal_aero_wateruptake


